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The nature of the driven steady states of many inter- 
acting particles, and the transitions between them, is a 
topic of much active interest. As with equilibrium sys- 
tems, the use of lattice models, in which the degrees of 
freedom sit on the sites of a discrete grid, has led to ana- 
lytical simplifications and greater accuracy in numerical 
simulations, as compared to contimmm models, [ij. Here 
we consider lattice models applied to driven two dimen- 
sional (2D) charges on a triangular grid, as a model for 
vortices in a honeycomb superconducting network. We 
use two distinct lattice gas dynamics, both intended to 
model the overdamped diffusive limit: (i) the commonly 
used driven diffusive Metropolis Monte Carlo (MC) 
and its modification to a (ii) driven diffusive continuous 
time Monte Carlo 0, 0| ■ We believe this is the first ap- 
plication of continuous time MC in the context of driven 
diffusive problems. We find that the steady states are 
qualitatively different for the two dynamics, and that the 
latter gives the more physically reasonable results. We 
find that finite size effects can be subtle at low tempera- 
ture. 

Our model is given by the Hamiltonian Q , 
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where rii = 0, 1 is the charge on site of a periodic 
triangular grid, — / is a fixed uniform background charge, 
and V{r) is the 2D lattice Coulomb interaction as defined 
in Ref.[5j for a triangular grid with periodic boundary 
conditions. We take as grid basis vectors fii — x, and 
(12 = '^x + the grid size is Li in direction a^, and 
the grid sites are r = miSi -I- m2Ci2, mi = 0, . . . , Li — 1. 
Neutrality requires a fixed number of charges, Xi "-i = 
Nc = /L1L2. In this work we use a charge density of 
/ = 1/25. In equilibrium, this model is characterized by 
a single first order melting transition at Tm — 0.009 from 
a triangular commensurate pinned solid with long range 
translational order to an ordinary liquid . 

We now consider behavior in a uniform driving force, 
F = Fx, parallel to the grid direction di. We con- 



sider two different dynamics, both involving single parti- 
cle moves only. 

(i) Driven diffusive Metropohs Monte Carlo (DDMC) 
0: At each step of the simulation, a charge ni = 1 
is selected at random, and moved a distance Ar to a 
nearest neighbor site. If 03 = ai — 0,2, then Ar is chosen 
randomly from the six possibilities, ±0^, i — 1,2,3. If 
Tioid and Tincw give the interaction energy before and 
after the move, one computes, 



AL; = T^ncw - Wold - F • Ar 
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where the last term is the work done by the force on the 
moving charge. One accepts or rejects this move accord- 
ing to the usual Metropolis MC algorithm. One pass of 
Nc steps equals one unit of simulation time. Statisti- 
cal averages are computed averaging over the generated 
configurations as in ordinary MC. 

(ii) Driven diffusive continuous time Monte Carlo 
(CTMC) HQ: At each step of the simulation, one con- 
siders the possible move of each charge Ui = 1 in each of 
the six directions, a = ±0^, computing the energy change 
AL'iQ, of each move according to Eq.lj^J. We take the rate 
for a particular move ia to be. 
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where I/Wq sets the unit of time. The total rate for all 
single particle moves is then Wtot = Xia ^ia- We decide 
which move to make by sampling the probability distri- 
bution Pia = Wia/Wtot, and then update the simulation 
clock by At = 1/Wtot- Averages of an observable O are 
computed as. 



{O) 
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where s labels the steps of the simulation, Og is the value 
of O in the configuration at time ts, Ai^ = ij+i ~ts, and 
T = Xs is the total time of the simulation. 

The CTMC, originally introduced as the "n-fold way" 
for spin models j3| , owes its name to the continuous vari- 
ations in the time steps At^, which vary as the configu- 
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ration changes throughout the simulation. It is a rejec- 
tionless algorithm designed to speed up excitation over 
energy barriers at low temperatures; rather than waste 
many rejected moves until a rare acceptance takes one 
up an energy barrier, the energy barriers AE themselves 
set the time scale for each move, which then happens in a 
single step. Simulation clock times can vary over orders 
of magnitude as T varies. 

In CTMC, there are many possible choices for the rates 
that will obey local detailed balance. It can be shown 
that the rates of Eq. Q lead to ordinary Langevin dynam- 
ics in the limit AEi^/T << 1. Our simulations, however, 
are generally in the opposite limit AEia/T > 1. To see 
what physical limit this corresponds to, consider a single 
particle on a one dimensional grid in a driving force F. 
A simple calculation gives for the average velocity of the 
particle, (v) = I^tot tanh(F/2r) = 2Wo sinh(F/2T). If 
we interpret the grid sites as the minima of a periodic 
pinning potential U(r) in the continuum, with energy 
barrier Uq, then Wq ^ e~^°/'^ and the above velocity 
then describes motion in such a periodic potential in the 
limit F <Uq 0. CTMC thus appears to describe the 
limit where motion is due to thermal activation over bar- 
riers; it is unclear if it can describe the very large drive 
limit, where the washboard potential C/(r) — F • r loses 
its local minima parallel to F. This large drive limit has 
been the subject of numerous theoretical d, 0, 0| and 
numerical |llL Il2l Il3l [iJI works for the case of random 
pinning. 

In this paper we consider just the structural properties 
of the steady state. These are given by the structure 
function, 
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where (jk — X^i e*'' '"' (n^ — /) is the Fourier transform of 
the charge distribution, k — fcibi 4- fc2b2 are the allowed 
wave vectors, with bi = 2'kx — and b2 — the 
basis vectors of the reciprocal lattice to the grid, and 
ki = li/Li with Zi = 0, . . . ii-i. We also consider the real 
space correlations. 



C(mi, fe) ^Y[Y. e-"''"''S{ki,k2) 
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as well as C{ki, TO2) defined similarly, and the 6-fold (hex- 
atic) orientational order parameter ($6)j where 
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In the above, the first sum is over all charges rii = 1, the 
second sum is over all charges nj = 1 that are nearest 
neighbors of Ui (as determined by a Delaunay triangu- 
lation), Zi is the number of these nearest neighbors, and 
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FIG. 1: (a) Plot of structure function transverse peak height 
^(Ki), Ki = |b2, vs. T, and (b) hexatic order parameter 
($6) vs. T, at fixed F = 0.1s, for both DDMC and CTMC 
algorithms. 



9ij is the angle of the bond from rii to Uj with respect to 
the X axis. 

The direct computation of S(k) and (^g) as in Eq.l^J 
is too costly as it involves lengthy calculations at each 
step of the simulation. Instead we approximate the in- 
tegral in Eq.Q by a Monte Carlo evaluation, choosing 
-^config — 1000 random times uniformly distributed over 
the simulation interval t G [0, r] and averaging over the 
configurations at these times only. 

We now present our results. Starting in the ground 
state at T = 0, the charge lattice remains pinned until 
the driving force F exceeds the change in interaction en- 
ergy associated with moving one charge forward one grid 
space parallel to F. This defines the T = critical force, 
Fc = 0.063, for both DDMC and CTMC. Our simula- 
tions are carried out for fixed F — 0.1 as we vary T. Our 
results reported here are for a system size of Li = 500 
and L2 = 60, with Nc = 1200 charges. The reason for 
such an extreme aspect ratio will be discussed later. At 
F = 0, T < Tm, the system forms an ordered triangular 
charge solid, and S'(k) has sharp Bragg peaks at recip- 
rocal lattice vectors {K}. Let Ki = b2/5 be smallest 
non-zero reciprocal lattice vector directed transverse to 
F. In Fig.IBi we plot S'(Ki) vs. T, at = 0.1, for both 
DDMC and CTMC. In Fig.^ we plot (^g) vs. T. The 
results for DDMC show no structure whatever for the 
moving steady state. A plot in Fig.|2t of the full S'(k) at 
T = 0.003 shows an isotropic liquid. In fact, we find with 
DDMC that once the charged solid depins from the grid, 
the moving steady state is an isotropic liquid virtually 
everywhere in the F — T plane. 

To see why this is so, consider a large F » Fc at 
T = 0, starting from an ordered solid. The DDMC picks 
a charge at random, then picks a direction to move it 
in at random; the move is accepted only if it lowers the 
energy, i.e. if the charge advances in the direction of F. 
Since only 3 of the 6 possible directions do so, the move 
is accepted with probability 1/2. If the work done by the 
force dominates the interaction energy (as it should for 



3 




FIG. 2: Structure function S{k), for k in the first Brillouin 
zone, at force F = 0.14 for (a) T = 0.003 with DDMC algo- 
rithm; and (b) T = 0.003, (c) T = 0.004, (d) T = 0.007 with 
CTMC algorithm. 



F >> Fc) then after one pass a random half of all charges 
have moved forward. On the next pass, a different ran- 
dom half move forward. After several passes, the charges 
are completely disordered. Although this argument as- 
sumed F >> Fc, we find that at T = the charges melt 
to a liquid at all F > Fc- The randomness of choosing 
proposed moves thus has a dramatic effect on the steady 
state order. In contrast, in CTMC, moves are chosen 
according to a probabilistic distribution which sharpens 
dramatically as T decreases. Unlikely moves are almost 
never chosen, and favorable moves are almost always cho- 
sen. The result, shown in Fig.^ is more structure in the 
moving steady states. The rest of this work, therefore, 
focuses on CTMC. 

For CTMC, the results of Fig. ^ show a discontinuous 
melting transition at Tm = 0.00325. In Fig.l^b— d we 
show representative plots of S'(k) above and below Tm. 
We first consider T > T^. Although the plot of S'(k) 
at T = 0.004 in Fig.|3; shows sharp peaks at the re- 
ciprocal lattice vectors of the ordered charge lattice,the 
magnitude of these peaks is greatly reduced from those 
at T < Tin (see Fig.|5j3). We have carried out simula- 
tions for a larger system, L2 = 120, Li = 500 and found 
S'(k) to remain unchanged. This lack of finite size de- 
pendence in S'(k) indicates that the system is a liquid 
with short ranged translational order. The peaks in S{k) 
result from large but finite correlation lengths. Similar 
behavior was seen in simulations of vortices in a square 
Josephson junction array with random pinning the 
prominent peaks at the transverse wavevectors along the 
ky axis led those authors to denote this state as a "short 




ranged smectic" 



FIG. 3: (a) Transverse real space correlation C{ki — 
0,7712) vs. m2, and (b) longitundinal real space corre- 
lation C{m\,k2 — 0) vs. mi. Solid lines are fits to 
g-(i-™.)/€ g^jjjj determine the correlation lengths 
and Curves from top to bottom are for T = 

0.00325, 0.0035, 0.004, 0.0045, 0.005, 0.0055, 0.006, 0.007; F = 
0.1. Representative values for ^± and ^|| are shown. 



To estimate the correlations lengths we plot in Fig.|31i,b 
C(fci — 0, vs. 7712 and 0(1711, k2 = 0) vs. mi; the 
first gives the decay of correlations in real space along 02 
(averaged over the direction ai), while the second gives 
the decay of correlations along ai, parallel to the force 
F (averaged over 02). We plot only the values at integer 
multiples of the average particle spacing a„ = l/-// = 5; 
these define the upper envelop of the damped oscillat- 
ing correlations. Fitting the data to the simple periodic 
decay C ~ e"'"/'^ -I- e"^^""'^/^ (where we use Li or L2 
as appropriate) gives the correlation lengths perpendic- 
ular, ^j^, and parallel, ^||, to the driving force F. Repre- 
sentative values for and ^|| are shown in Fig. |31 with 
^ 2(^11 near Tm. For Fig.l^t) our fit is only to points 

Although the liquid above Tm lacks translational or- 
der, Fig.^ shows that hexatic orientational order grows 
for T < T* ^ 0.007. Similar hexatic liquids have been 
reported in continuum simulations with random pinning 
|12l |. In our case, the periodic triangular grid breaks ro- 
tational symmetry and in principle gives finite hexatic 
order at any T. It is unclear whether the onset of grow- 
ing hexatic order at T* is just this grid induced effect, 
increasing as the correlations lengths grow larger than the 
interparticle spacing, ^ > a« , or whether it is a crossover 
remnant of what might be a true hexatic transition in 
another geometry. 

Next we consider T < T^. S{k.) for T = 0.003 < T„i is 
shown in Fig.lJt. First we note that the peaks at K on 
the ky axis (at fci = 0, fc2 = 1/5, 2/5) are sharp (5-function 
Bragg peaks of height S'(K) ~ iVc = 1200. We have 
computed the height of these peaks for smaller size sys- 
tems and find that they scale ~ Nc- This indicates that 
this state has long range smectic order; the particles are 
confined to periodically spaced channels parallel to the 
driving force. Next, we note that the peaks at finite k^ 
(at ki = 1/5,2/5) are essentially (5-functions in k^. Sim- 
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FIG. 4: (a) Transverse correlation between the smectic chan- 
nels, C{ki — 1/5, 7712) vs. m2, and (b) longitudinal cor- 
relation C(mi,fc2 — 0) vs. mi, in the smectic state for 
T = 0.003,0.00325 < T^, J' = 0.1. Solid lines are fits to 
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ulations for smaller size systems show that the height of 
these peaks scale ^ Li. This indicates that the parti- 
cles are perfectly ordered within each smectic channel. 
The finite width of these peaks, as ky varies, indicates 
that the ordered channels are randomly displaced with 
respect to each other, with a finite correlation length 
To determine this correlation length we plot in Fig.^ 
C(fci — 1/5,7712) vs. 7772 and fit to a periodic decay as 
earlier. Finally we return to the issue of the longitu- 
dinal order in the smectic state. Since the transverse 
order among the smectic channels is short ranged, with 
finite we can regard the channels as decoupled one 
dimensional systems. Thus only short ranged longitu- 
dinal order should be expected. We have investigated 
this issue by carrying out detailed finite size analyses on 
smaller size systems. We conclude that the smectic does 
in fact have a finite longitudinal correlation length , 
but that ^ Li; we find that when the system length 
is increased so that < Li, the smectic becomes un- 
stable to the liquid. In Fig.^ we plot the longitudinal 
correlation C{mi,k2 = 0) vs. mi for T = 0.003,0.00325 
just below Tm. Fitting to the periodic exponential de- 
cay assumed earlier, we find ~ 500. For a smaller 
length, Li — 120, the smectic persisted up to the higher 
T = 0.006. Our desire to supress to low tempera- 
tures, so as to see growing correlations in the liquid, was 
the reason we chose Li = 500 for the simulations re- 
ported on here. While the smectic thus disappears in the 
true thermodynamic limit, since grows exponentially 
as T decreases, the smectic will ultimately appear in a 
finite sized system at sufficiently low T. We find that 
once > Li, the smectic is the stable state of the sys- 
tem; for Li = L2 = 120, we have suceeded in cooling into 
the smectic from the disordered liquid. Our observation 
of a smectic state on finite length scales, which becomes 
unstable to a liquid on large length scales, agrees in part 
with arguments in Ref. 



Interacting 2D vortices in a periodic potential at finite 
temperature have been simulated by several others using 
continuum dynamics. The molecular dynamic simula- 
tions of Reichhardt and Zimanyi and Carneiro 
used square periodic pins embedded in a flat continuum, 
with a number of vortices equal to, or greater than, the 
number of pins. Such models cannot be well described 
by our discrete grid. Closer to our model is that of Mar- 
coni and Domfnguez [l^ who simulate the RSJ dynam- 
ics of a vortex density / — 1/25 in a square Josephson 
array, and find an ordered moving vortex lattice. How- 
ever in their case, the energy to move a single vortex 
forward from its ground state position is AEi ~ 0.34, 
whereas the energy barrier between cells of the array is 
Uq ~ 0.12. The parameters of our simulations, which as- 
sume AEi < F < Ua (see discussion preceding Eq. Q), 
are therefore in a more strongly pinned limit outside the 
range of their model. It therefore remains for future in- 
vestigation to test if the results of the CTMC method 
agree with that of continuum models in the correspond- 
ing limit. 

We wish to thank D. Domfnguez, M. C. Marchetti and 
A. A. Middleton for helpful discussions. This work was 
supported by DOE grant DE-FG02-89ER14017 and NSF 
grant PHY-9987413. 
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